//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#include <buola/image/cpolygon.h>
#include <queue>

namespace buola { namespace img {

/* Compare the slope of a to the slope of b, returning 1, 0, -1 if the
 * slope a is respectively greater than, equal to, or less than the
 * slope of b.
 *
 * For each edge, consider the direction vector formed from:
 *
 *        top -> bottom
 *
 * which is:
 *
 *        (dx, dy) = (line.p2.x - line.p1.x, line.p2.y - line.p1.y)
 *
 * We then define the slope of each edge as dx/dy, (which is the
 * inverse of the slope typically used in math instruction). We never
 * compute a slope directly as the value approaches infinity, but we
 * can derive a slope comparison without division as follows, (where
 * the ? represents our compare operator).
 *
 * 1.           slope(a) ? slope(b)
 * 2.            adx/ady ? bdx/bdy
 * 3.        (adx * bdy) ? (bdx * ady)
 *
 * Note that from step 2 to step 3 there is no change needed in the
 * sign of the result since both ady and bdy are guaranteed to be
 * greater than or equal to 0.
 *
 * When using this slope comparison to sort edges, some care is needed
 * when interpreting the results. Since the slope compare operates on
 * distance vectors from top to bottom it gives a correct left to
 * right sort for edges that have a common top point, (such as two
 * edges with start events at the same location). On the other hand,
 * the sense of the result will be exactly reversed for two edges that
 * have a common stop point.
 */
static inline bool slope_less(const CPolygon::SEdge &p1,const CPolygon::SEdge &p2)
{
    /* XXX: We're assuming here that dx and dy will still fit in 32
     * bits. That's not true in general as there could be overflow. We
     * should prevent that before the tessellation algorithm
     * begins.
     */

    auto lDX1=p1.mLine.m1.x-p1.mLine.m2.x;
    auto lDX2=p2.mLine.m1.x-p2.mLine.m2.x;

    /* Since the dy's are all positive by construction we can fast
     * path several common cases.
     */

    /* First check for vertical lines. */
    if (!lDX1)
        return lDX2>nZero;
    if (!lDX2)
        return lDX1<nZero;

    /* Then where the two edges point in different directions wrt x. */
    if (lDX1>nZero)
    {
        if(lDX2<nZero) return false;
    }
    else if(lDX2>nZero) return true;

    /* Finally we actually need to do the general comparison. */
    auto lDY1=p1.mLine.m1.y-p1.mLine.m2.y;
    auto lDY2=p2.mLine.m1.y-p2.mLine.m2.y;

    return lDX1*lDY2<lDX2*lDY1;
}

/* Compute the intersection of two edges. The result is provided as a
 * coordinate pair of 128-bit integers.
 *
 * Returns %CAIRO_BO_STATUS_INTERSECTION if there is an intersection
 * that is within both edges, %CAIRO_BO_STATUS_NO_INTERSECTION if the
 * intersection of the lines defined by the edges occurs outside of
 * one or both edges, and %CAIRO_BO_STATUS_PARALLEL if the two edges
 * are exactly parallel.
 *
 * Note that when determining if a candidate intersection is "inside"
 * an edge, we consider both the infinitesimal shortening and the
 * infinitesimal tilt rules described by John Hobby. Specifically, if
 * the intersection is exactly the same as an edge point, it is
 * effectively outside (no intersection is returned). Also, if the
 * intersection point has the same
 */
static inline CPoint_x edge_intersection(const CPolygon::SEdge &p1,const CPolygon::SEdge &p2)
{
    auto lDX1=p1.mLine.m1.x-p1.mLine.m2.x;
    auto lDX2=p2.mLine.m1.x-p2.mLine.m2.x;
    auto lDY1=p1.mLine.m1.y-p1.mLine.m2.y;
    auto lDY2=p2.mLine.m1.y-p2.mLine.m2.y;


    auto lDenDet=lDX1*lDY2-lDY1*lDX2;

     /* Q: Can we determine that the lines do not intersect (within range)
      * much more cheaply than computing the intersection point i.e. by
      * avoiding the division?
      *
      *   X = ax + t * adx = bx + s * bdx;
      *   Y = ay + t * ady = by + s * bdy;
      *   ∴ t * (ady*bdx - bdy*adx) = bdx * (by - ay) + bdy * (ax - bx)
      *   => t * L = R
      *
      * Therefore we can reject any intersection (under the criteria for
      * valid intersection events) if:
      *   L^R < 0 => t < 0, or
      *   L<R => t > 1
      *
      * (where top/bottom must at least extend to the line endpoints).
      *
      * A similar substitution can be performed for s, yielding:
      *   s * (ady*bdx - bdy*adx) = ady * (ax - bx) - adx * (ay - by)
      */
    auto R = (p2.mLine.m1.y-p1.mLine.m1.y)*lDX2-(p2.mLine.m1.x-p1.mLine.m1.x)*lDY2;
    if(lDenDet<nZero)
    {
        if(lDenDet>=R) return nUndefined;
    }
    else
    {
        if(lDenDet<=R) return nUndefined;
    }

    R = (p1.mLine.m1.y-p2.mLine.m1.y)*lDX1-(p1.mLine.m1.x-p2.mLine.m1.x)*lDY1;
    if(lDenDet<nZero)
    {
        if(lDenDet>=R) return nUndefined;
    }
    else
    {
        if(lDenDet<=R) return nUndefined;
    }

    auto lDet1=p1.mLine.m1.x*p1.mLine.m2.y-p1.mLine.m2.x*p1.mLine.m1.y;
    auto lDet2=p2.mLine.m1.x*p2.mLine.m2.y-p2.mLine.m2.x*p2.mLine.m1.y;

    ///\todo ensure that this is numerically accurate by performing it with 128 bits and doing much more
    ///\todo we should also check and store whether it is exact
    CPoint_x lIntersection={(lDet1*lDX2-lDX1*lDet2)/lDenDet,(lDet1*lDY2-lDY1*lDet2)/lDenDet};

    if(!p1.mLine.ContainsIntersection(lIntersection)) return nUndefined;
    if(!p2.mLine.ContainsIntersection(lIntersection)) return nUndefined;

    return lIntersection;
}

bool slope_equal(const CPolygon::SEdge &p1,const CPolygon::SEdge &p2)
{
    /* XXX: We're assuming here that dx and dy will still fit in 32
     * bits. That's not true in general as there could be overflow. We
     * should prevent that before the tessellation algorithm
     * begins.
     */
    auto lDX1=p1.mLine.m2.x-p1.mLine.m1.x;
    auto lDX2=p2.mLine.m2.x-p2.mLine.m1.x;
    auto lDY1=p1.mLine.m1.y-p1.mLine.m2.y;
    auto lDY2=p2.mLine.m1.y-p2.mLine.m2.y;

    return lDX1*lDY2==lDX2*lDY1;
}

bool edges_colinear(const CPolygon::SEdge &p1,const CPolygon::SEdge &p2)
{
    bool lX1=p1.mLine.m1.x==p2.mLine.m1.x;
    bool lX2=p1.mLine.m2.x==p2.mLine.m2.x;
    bool lY1=p1.mLine.m1.y==p2.mLine.m1.y;
    bool lY2=p1.mLine.m2.y==p2.mLine.m2.y;

    if(lX1&&lX2&&lY1&&lY2) return true;

    if(!slope_equal(p1,p2)) return false;

    if(lX1||lX2||lY1||lY2)
        return (lX1&&lY1)||(lX2&&lY2);

    return p1.mLine.m1.x==p2.mLine.GetXForY(p1.mLine.m1.y);
}

/* Adds the trapezoid, if any, of the left edge to the #cairo_traps_t */
void edge_end_trap(CPolygon::SEdge *pLeft,const CFixed<> &pBottom,std::vector<CTrap> &pTraps)
{
    CPolygon::SDeferredTrap *lTrap=&pLeft->mDeferredTrap;

    /* Only emit (trivial) non-degenerate trapezoids with positive height. */
    if(lTrap->mTop<pBottom)
    {
        CTrap lNew;
        lNew.mLeft=pLeft->mLine;
        lNew.mRight=lTrap->mRight->mLine;
        lNew.mTop=lTrap->mTop;
        lNew.mBottom=pBottom;

        pTraps.push_back(lNew);
    }

    lTrap->mRight=nullptr;
}

bool operator<(const CPolygon::SEvent &pLH,const CPolygon::SEvent &pRH)
{
    if(pLH.mPoint.y!=pRH.mPoint.y)
        return pLH.mPoint.y<pRH.mPoint.y;
    if(pLH.mPoint.x!=pRH.mPoint.x)
        return pLH.mPoint.x<pRH.mPoint.x;

    if(pLH.mType!=pRH.mType)
        return pLH.mType<pRH.mType;

    return &pLH<&pRH;
}

struct CBOEventGreater
{
    bool operator()(const CPolygon::SEvent *p1,const CPolygon::SEvent *p2)
    {
        return *p2<*p1;
    }
};

struct CBOEventLess
{
    bool operator()(const CPolygon::SEvent *p1,const CPolygon::SEvent *p2)
    {
        return *p1<*p2;
    }
};

class CPolygon::SSweepLine
{
public:
    SSweepLine()
        :   mHead(nullptr)
        ,   mStopped(nullptr)
        ,   mCurrentY(std::numeric_limits<CFixed<>>::min())
        ,   mCurrentEdge(nullptr)
    {}

    int CompareEdges(const SEdge &p1,const SEdge &p2)
    {
        ///\todo this
        using std::max;
        using std::min;
        
        if(p1.mLine!=p2.mLine)
        {
            if(max(p1.mLine.m1.x,p1.mLine.m2.x)<min(p2.mLine.m1.x,p2.mLine.m2.x))
                return -1;

            if(min(p1.mLine.m1.x,p1.mLine.m2.x)>max(p2.mLine.m1.x,p2.mLine.m2.x))
                return 1;

            auto lX1=p1.mLine.GetXForY(mCurrentY);
            auto lX2=p2.mLine.GetXForY(mCurrentY);

            if(lX1<lX2)
                return -1;
            if(lX2<lX1)
                return 1;

            ///\todo slope comparison
        }

        /* We've got two collinear edges now. */
        return (p2.mLine.mBottom-p1.mLine.mBottom).Get();
    }

    void Insert(SEdge *pEdge)
    {
        if(mCurrentEdge)
        {
            SEdge *lPrev,*lNext;

            int lCmp=CompareEdges(*mCurrentEdge,*pEdge);

            if(lCmp<0)
            {
                lPrev=mCurrentEdge;
                lNext=lPrev->mNext;
                while(lNext&&CompareEdges(*lNext,*pEdge)<0)
                {
                    lPrev=lNext;
                    lNext=lPrev->mNext;
                }

                lPrev->mNext=pEdge;
                pEdge->mPrev=lPrev;
                pEdge->mNext=lNext;
                if(lNext)
                    lNext->mPrev=pEdge;
            }
            else if(lCmp>0)
            {
                lNext=mCurrentEdge;
                lPrev=lNext->mPrev;
                while(lPrev&&CompareEdges(*lPrev,*pEdge)>0)
                {
                    lNext=lPrev;
                    lPrev=lNext->mPrev;
                }

                lNext->mPrev=pEdge;
                pEdge->mNext=lNext;
                pEdge->mPrev=lPrev;
                if(lPrev)
                    lPrev->mNext=pEdge;
                else
                    mHead=pEdge;
            }
            else
            {
                lPrev=mCurrentEdge;
                pEdge->mPrev=lPrev;
                pEdge->mNext=lPrev->mNext;
                if(lPrev->mNext)
                    lPrev->mNext->mPrev=pEdge;
                lPrev->mNext=pEdge;
            }
        }
        else
        {
            mHead=pEdge;
            pEdge->mNext=nullptr;
            pEdge->mPrev=nullptr;
        }

        mCurrentEdge=pEdge;
    }

    void Delete(SEdge *pEdge)
    {
        if(pEdge->mPrev)
            pEdge->mPrev->mNext=pEdge->mNext;
        else
            mHead=pEdge->mNext;

        if(pEdge->mNext)
            pEdge->mNext->mPrev=pEdge->mPrev;

        if(mCurrentEdge==pEdge)
            mCurrentEdge=pEdge->mPrev?pEdge->mPrev:pEdge->mNext;
    }

    void Swap(SEdge *pLeft,SEdge *pRight)
    {
        if(pLeft->mPrev)
            pLeft->mPrev->mNext=pRight;
        else
            mHead=pRight;

        if(pRight->mNext)
            pRight->mNext->mPrev=pLeft;

        pLeft->mNext=pRight->mNext;
        pRight->mNext=pLeft;
        pRight->mPrev=pLeft->mPrev;
        pLeft->mPrev=pRight;
    }

    SEdge *mHead;
    SEdge *mStopped;
    CFixed<> mCurrentY;
    SEdge *mCurrentEdge;
};

class CPolygon::SEventQueue
{
public:
    SEventQueue(std::vector<SEvent> &pStart)
        :   mStart(pStart.size())
    {
        for(int i=0;i<pStart.size();i++)
        {
            mStart[i]=&pStart[i];
        }

        std::sort(mStart.begin(),mStart.end(),CBOEventLess());

        mStartPos=0;
    }
    
    SEvent *Dequeue()
    {
        if(mQueue.empty())
        {
            if(mStartPos<mStart.size())
                return mStart[mStartPos++];
            return nullptr;
        }
        
        SEvent *lEvent=mQueue.top();

        if(mStartPos<mStart.size())
        {
            SEvent *lCmp=mStart[mStartPos];
        
            if (lCmp&&*lCmp<*lEvent)
            {
                mStartPos++;
                return lCmp;
            }
        }

        mQueue.pop();
        return lEvent;
    }

    void InsertStop(SEdge *pEdge)
    {
        SEvent *lEvent=new SEvent;
        lEvent->mType=BO_STOP;
        lEvent->mPoint.y=pEdge->mLine.mBottom;
        lEvent->mPoint.x=pEdge->mLine.GetXForY(lEvent->mPoint.y);
        lEvent->mEdge.mNext=pEdge;

        mQueue.push(lEvent);
    }

    void InsertIfIntersectBelowCurrentY(SEdge *pLeft,SEdge *pRight)
    {
        using std::max;
        using std::min;

        if(max(pLeft->mLine.m1.x,pLeft->mLine.m2.x)<=min(pRight->mLine.m1.x,pRight->mLine.m2.x))
            return;

        if(pLeft->mLine==pRight->mLine)
            return;

        /* The names "left" and "right" here are correct descriptions of
        * the order of the two edges within the active edge list. So if a
        * slope comparison also puts left less than right, then we know
        * that the intersection of these two segments has already
        * occurred before the current sweep line position. */
        if (!slope_less(*pRight,*pLeft))
            return;

        CPoint_x lIntersection=edge_intersection(*pLeft,*pRight);

        if(lIntersection.Undefined()) return;

        SEvent *lEvent=new SEvent;
        lEvent->mType=BO_INTERSECTION;
        lEvent->mPoint=lIntersection;
        lEvent->mEdge.mNext=pLeft;
        lEvent->mEdge.mPrev=pRight;

        mQueue.push(lEvent);
    }

private:
    std::priority_queue<SEvent*,std::vector<SEvent*>,CBOEventGreater> mQueue;
    std::vector<SEvent*> mStart;
    int mStartPos;
};

/* Start a new trapezoid at the given top y coordinate, whose edges
 * are `edge' and `edge->next'. If `edge' already has a trapezoid,
 * then either add it to the traps in `traps', if the trapezoid's
 * right edge differs from `edge->next', or do nothing if the new
 * trapezoid would be a continuation of the existing one. */
static inline void start_or_continue_trap(CPolygon::SEdge *pLeft,CPolygon::SEdge *pRight,const CFixed<> &pTop,
                                          std::vector<CTrap> &pTraps)
{
    if(pLeft->mDeferredTrap.mRight==pRight)
        return;

    if(pLeft->mDeferredTrap.mRight)
    {
        if(edges_colinear(*pLeft->mDeferredTrap.mRight,*pRight))
        {
            pLeft->mDeferredTrap.mRight=pRight;
            return;
        }

        edge_end_trap(pLeft,pTop,pTraps);
        return;
    }

    if(!edges_colinear(*pLeft,*pRight))
    {
        pLeft->mDeferredTrap.mTop=pTop;
        pLeft->mDeferredTrap.mRight=pRight;
    }
}

void active_edges_to_traps(CPolygon::SEdge *pPos,const CFixed<> &pTop,uint32_t pWinding,std::vector<CTrap> &pTraps)
{
    CPolygon::SEdge *lLeft=pPos;
    int lInOut=0;

    while(pPos)
    {
        if(pPos!=lLeft&&pPos->mDeferredTrap.mRight)
        {
            /* XXX It shouldn't be possible to here with 2 deferred traps
             * on colinear edges... See bug-bo-rictoz.
             */
            if(!lLeft->mDeferredTrap.mRight&&edges_colinear(*lLeft,*pPos))
            {
                /* continuation on left */
                lLeft->mDeferredTrap=pPos->mDeferredTrap;
                pPos->mDeferredTrap.mRight=nullptr;
            }
            else
            {
                edge_end_trap(pPos,pTop,pTraps);
            }
        }

        lInOut+=pPos->mDirection;
        if((lInOut&pWinding)==0)
        {
            /* skip co-linear edges */
            if(!pPos->mNext||!edges_colinear(*pPos,*pPos->mNext))
            {
                start_or_continue_trap(lLeft,pPos,pTop,pTraps);

                lLeft = pPos->mNext;
            }
        }

        pPos = pPos->mNext;
    }
}

std::vector<CTrap> CPolygon::TesselateToTraps()
{
    std::vector<CTrap> lTraps;

    uint32_t lFillRule=(uint32_t)-1;

    CPolygon::SEventQueue lQueue(mEdges);
    CPolygon::SSweepLine lSweepLine;

    while(CPolygon::SEvent *lEvent=lQueue.Dequeue())
    {
        if(lEvent->mPoint.y!=lSweepLine.mCurrentY)
        {
            for(CPolygon::SEdge *lE1=lSweepLine.mStopped;lE1;lE1=lE1->mNext)
            {
                if(lE1->mDeferredTrap.mRight)
                    edge_end_trap(lE1,lE1->mLine.mBottom,lTraps);
            }
            lSweepLine.mStopped=nullptr;

            active_edges_to_traps(lSweepLine.mHead,lSweepLine.mCurrentY,lFillRule,lTraps);

            lSweepLine.mCurrentY=lEvent->mPoint.y;
        }

        switch(lEvent->mType)
        {
        case CPolygon::BO_START:
            {
                CPolygon::SEdge *lE1=&lEvent->mEdge;

                lSweepLine.Insert(lE1);
                lQueue.InsertStop(lE1);

                /* check to see if this is a continuation of a stopped edge */
                /* XXX change to an infinitesimal lengthening rule */
                for(CPolygon::SEdge *lLeft=lSweepLine.mStopped;lLeft;lLeft=lLeft->mNext)
                {
                    if(lE1->mLine.mTop<=lLeft->mLine.mBottom&&edges_colinear(*lE1,*lLeft))
                    {
                        lE1->mDeferredTrap=lLeft->mDeferredTrap;
                        if(lLeft->mPrev)
                            lLeft->mPrev=lLeft->mNext; ///this doesn't seem right!!!
                        else
                            lSweepLine.mStopped=lLeft->mNext;
                        if(lLeft->mNext)
                            lLeft->mNext->mPrev=lLeft->mPrev;
                        break;
                    }
                }

                CPolygon::SEdge *lLeft=lE1->mPrev;
                CPolygon::SEdge *lRight=lE1->mNext;

                if(lLeft)
                    lQueue.InsertIfIntersectBelowCurrentY(lLeft,lE1);

                if(lRight)
                    lQueue.InsertIfIntersectBelowCurrentY(lE1,lRight);

                break;
            }

        case CPolygon::BO_STOP:
            {
                CPolygon::SEdge *lE1=lEvent->mEdge.mNext;
                delete lEvent; ///why? I guess it means really delete as in C++, because it was created

                CPolygon::SEdge *lLeft=lE1->mPrev;
                CPolygon::SEdge *lRight=lE1->mNext;

                lSweepLine.Delete(lE1);

                if(lE1->mDeferredTrap.mRight)
                {
                    lE1->mNext=lSweepLine.mStopped;
                    if(lSweepLine.mStopped)
                        lSweepLine.mStopped->mPrev=lE1;
                    lSweepLine.mStopped=lE1;
                    lE1->mPrev=nullptr;
                }

                if(lLeft&&lRight)
                    lQueue.InsertIfIntersectBelowCurrentY(lLeft,lRight);

                break;
            }

        case CPolygon::BO_INTERSECTION:
            {
                CPolygon::SEdge *lE1=lEvent->mEdge.mNext;
                CPolygon::SEdge *lE2=lEvent->mEdge.mPrev;

                delete lEvent;

                if(lE2!=lE1->mNext)
                    break;

                CPolygon::SEdge *lLeft=lE1->mPrev;
                CPolygon::SEdge *lRight=lE2->mNext;

                lSweepLine.Swap(lE1,lE2);

                if(lLeft)
                    lQueue.InsertIfIntersectBelowCurrentY(lLeft,lE2);
                if(lRight)
                    lQueue.InsertIfIntersectBelowCurrentY(lE1,lRight);
                break;
            }
        }
    }

    for(CPolygon::SEdge *lE1=lSweepLine.mStopped;lE1;lE1=lE1->mNext)
    {
        if(lE1->mDeferredTrap.mRight)
            edge_end_trap(lE1,lE1->mLine.mBottom,lTraps);
    }

    return lTraps;
}

/*namespace img*/ } /*namespace buola*/ }
